##########################################
# Replication file for "'Violence on Many Sides': Framing Effects on Protest and Support for Repression
# Edwards and Arnon, BJPS 2019
# Section 1: Israel Survey - Main article, then appendix (Figures 1&2 in Main, and Appendix Figure 2 and Tables 13-21)
#Section 2: US Survey - Main article, then appendix (Figures 3&4 in Main and Appendix Figure 3 and Tables 4-12)
#Section 3: Joint Data - Appendix Figure 1
########################################


#############################################################################
#END OF US SURVEY
#BEGINNING OF ISRAEL SURVEY
# MAKE SURE TO REMOVE ALL VARIABLES USING rm(list = ls()), 
# SINCE NAMES OF VARIABLE ARE IDENTICAL TO US SURVEY AT TIMES
############################################################################
rm(list = ls())
cat("\014")
set.seed(30030)


#1 load  data####
samp <- read.csv("Israel_Surveydata_final.csv")
#3 analysis####

#######################################################################
# Main hypothesis tests#
##Figures 1 and 2 in Main Text
############################################

#hypothesis 1: threat####

#difference in means, classification

ate.1 <- mean(samp$action[samp$threat == 1],na.rm=T) - mean(samp$action[samp$threat == 0],na.rm=T) #all
ate.2 <- mean(samp$action[samp$arm_conf == 1], na.rm=T) - mean(samp$action[samp$arm_enc == 1], na.rm = T) #within armed
ate.3 <- mean(samp$action[samp$un_conf == 1], na.rm=T) - mean(samp$action[samp$pr_nv == 1], na.rm = T) #within unarmed


#standard errors of ates

se.ate.1 <- sqrt(var(samp$action[samp$threat == 1],na.rm=T)/
                   nrow(subset(samp, samp$threat == 1 & !is.na(samp$action))) + 
                   var(samp$action[samp$threat == 0],na.rm=T)/
                   nrow(subset(samp, samp$threat == 0 & !is.na(samp$action))))

se.ate.2 <- sqrt(var(samp$action[samp$arm_conf == 1],na.rm=T)/
                   nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$action))) + 
                   var(samp$action[samp$arm_enc == 1],na.rm=T)/
                   nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$action))))

se.ate.3 <- sqrt(var(samp$action[samp$un_conf == 1],na.rm=T)/
                   nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$action))) + 
                   var(samp$action[samp$pr_nv == 1],na.rm=T)/
                   nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$action))))

#randomization inference

ri <- function(df, outcome, reps, n, m){
  
  library(randomizr)
  #create objects
  ATEs <- vector(mode="numeric")
  assns <- as.data.frame(matrix(NA,nrow=nrow(df),ncol=reps))
  dfs <- list()
  
  #generate random assignments
  for(i in 1:reps){
    assns[,i] <- complete_ra(N = n, m = m)}
  
  #remove duplicate assignments
  n=reps
  assns <- assns[!duplicated(lapply(assns, function(k) head(k, n)))]
  
  #calculate ATE for each assignment and store
  for(i in 1:ncol(assns)){dfs[[i]] <- cbind(outcome, assns[,i]) 
  ATEs[i] <- mean(dfs[[i]][,1][dfs[[i]][,2] == 1],na.rm=T) - mean(dfs[[i]][,1][dfs[[i]][,2] == 0],na.rm=T)}
  ATEs <- as.data.frame(sort(ATEs))
  
  #output
  return(ATEs)
}

seri.ate.1 <- ri(samp, samp$action, 1000, 1017, 467)
seri.ate.1.ci <- c(seri.ate.1[round(.025*nrow(seri.ate.1)),1], seri.ate.1[round(.975*nrow(seri.ate.1)), 1])

#difference in means, repressive response

ate.4 <- mean(samp$policy[samp$threat == 1],na.rm=T) - mean(samp$policy[samp$threat == 0],na.rm=T) #all
ate.5 <- mean(samp$policy[samp$arm_conf == 1], na.rm=T) - mean(samp$policy[samp$arm_enc == 1], na.rm = T) #within armed
ate.6 <- mean(samp$policy[samp$un_conf == 1], na.rm=T) - mean(samp$policy[samp$pr_nv == 1], na.rm = T) #within unarmed


#standard errors of ates

se.ate.4 <- sqrt(var(samp$policy[samp$threat == 1],na.rm=T)/
                   nrow(subset(samp, samp$threat == 1 & !is.na(samp$policy))) + 
                   var(samp$policy[samp$threat == 0],na.rm=T)/
                   nrow(subset(samp, samp$threat == 0 & !is.na(samp$policy))))

se.ate.5 <- sqrt(var(samp$policy[samp$arm_conf == 1],na.rm=T)/
                   nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$policy))) + 
                   var(samp$policy[samp$arm_enc == 1],na.rm=T)/
                   nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$policy))))

se.ate.6 <- sqrt(var(samp$policy[samp$un_conf == 1],na.rm=T)/
                   nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$policy))) + 
                   var(samp$policy[samp$pr_nv == 1],na.rm=T)/
                   nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$policy))))

#randomization inference

seri.ate.4 <- ri(samp, samp$policy, 1000, 1017, 467)
seri.ate.4.ci <- c(seri.ate.4[round(.025*nrow(seri.ate.4)),1], seri.ate.4[round(.975*nrow(seri.ate.4)), 1])


#regression

samp.armedonly <- subset(samp, samp$arm == 1)
samp.unarmedonly <- subset(samp, samp$arm == 0)

model.1 <- summary(lm(action ~ threat, data=samp))
model.2 <- summary(lm(action ~ threat, data=samp.armedonly))
model.3 <- summary(lm(action ~ threat, data=samp.unarmedonly))

model.4 <- summary(lm(policy ~ threat, data=samp))
model.5 <- summary(lm(policy ~ threat, data=samp.armedonly))
model.6 <- summary(lm(policy ~ threat, data=samp.unarmedonly))


#with covariates

model.1b <- summary(lm(action ~ threat + age + gender + income + education + religious_id + military, data=samp))
model.2b <- summary(lm(action ~ threat + age + gender + income + education + religious_id  + military, data=samp.armedonly))
model.3b <- summary(lm(action ~ threat + age + gender + income + education + religious_id  + military, data=samp.unarmedonly))

model.4b <- summary(lm(policy ~ threat + age + gender + income + education + religious_id  + military, data=samp))
model.5b <- summary(lm(policy ~ threat + age + gender + income + education + religious_id  + military, data=samp.armedonly))
model.6b <- summary(lm(policy ~ threat + age + gender + income + education + religious_id  + military, data=samp.unarmedonly))


#hypothesis 2: arms####

#difference in means, classification

ate.10 <- mean(samp$action[samp$arm == 1],na.rm=T) - mean(samp$action[samp$arm == 0],na.rm=T) #all
ate.11 <- mean(samp$action[samp$arm_conf == 1], na.rm=T) - mean(samp$action[samp$un_conf == 1], na.rm = T) #within threatening
ate.12 <- mean(samp$action[samp$arm_enc == 1], na.rm=T) - mean(samp$action[samp$pr_nv == 1], na.rm = T) #within unthreatening

#standard errors of ates

se.ate.10 <- sqrt(var(samp$action[samp$arm == 1],na.rm=T)/
                    nrow(subset(samp, samp$arm == 1 & !is.na(samp$action))) + 
                    var(samp$action[samp$arm == 0],na.rm=T)/
                    nrow(subset(samp, samp$arm == 0 & !is.na(samp$action))))

se.ate.11 <- sqrt(var(samp$action[samp$arm_conf == 1],na.rm=T)/
                    nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$action))) + 
                    var(samp$action[samp$un_conf == 1],na.rm=T)/
                    nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$action))))

se.ate.12 <- sqrt(var(samp$action[samp$arm_enc == 1],na.rm=T)/
                    nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$action))) + 
                    var(samp$action[samp$pr_nv == 1],na.rm=T)/
                    nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$action))))

#randomization inference

seri.ate.10 <- ri(samp, samp$action, 1000, 1017, 466)
seri.ate.10.ci <- c(seri.ate.10[round(.025*nrow(seri.ate.10)),1], seri.ate.10[round(.975*nrow(seri.ate.10)), 1])


#difference in means, repressive response

ate.13 <- mean(samp$policy[samp$arm == 1],na.rm=T) - mean(samp$policy[samp$arm == 0],na.rm=T) #all
ate.14 <- mean(samp$policy[samp$arm_conf == 1], na.rm=T) - mean(samp$policy[samp$un_conf == 1], na.rm = T) #within threatening
ate.15 <- mean(samp$policy[samp$arm_enc == 1], na.rm=T) - mean(samp$policy[samp$pr_nv == 1], na.rm = T) #within unthreatening


#standard errors of ates

se.ate.13 <- sqrt(var(samp$policy[samp$arm == 1],na.rm=T)/
                    nrow(subset(samp, samp$arm == 1 & !is.na(samp$policy))) + 
                    var(samp$policy[samp$arm == 0],na.rm=T)/
                    nrow(subset(samp, samp$arm == 0 & !is.na(samp$policy))))

se.ate.14 <- sqrt(var(samp$policy[samp$arm_conf == 1],na.rm=T)/
                    nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$policy))) + 
                    var(samp$policy[samp$un_conf == 1],na.rm=T)/
                    nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$policy))))

se.ate.15 <- sqrt(var(samp$policy[samp$arm_enc == 1],na.rm=T)/
                    nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$policy))) + 
                    var(samp$policy[samp$pr_nv == 1],na.rm=T)/
                    nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$policy))))

#randomization inference

seri.ate.13 <- ri(samp, samp$policy, 1000, 1017, 466)
seri.ate.13.ci <- c(seri.ate.13[round(.025*nrow(seri.ate.13)),1], seri.ate.13[round(.975*nrow(seri.ate.13)), 1])


#regression

samp.threatonly <- subset(samp, samp$threat == 1)
samp.unthreatonly <- subset(samp, samp$threat == 0)

model.10 <- summary(lm(action ~ arm, data=samp))
model.11 <- summary(lm(action ~ arm, data=samp.threatonly))
model.12 <- summary(lm(action ~ arm, data=samp.unthreatonly))

model.13 <- summary(lm(policy ~ arm, data=samp))
model.14 <- summary(lm(policy ~ arm, data=samp.threatonly))
model.15 <- summary(lm(policy ~ arm, data=samp.unthreatonly))

#with covariates

model.10b <- summary(lm(action ~ arm + age + gender + income + education + religious_id  + military, data=samp))
model.11b <- summary(lm(action ~ arm + age + gender + income + education + religious_id  + military, data=samp.threatonly))
model.12b <- summary(lm(action ~ arm + age + gender + income + education + religious_id  + military, data=samp.unthreatonly))

model.13b <- summary(lm(policy ~ arm + age + gender + income + education + religious_id  + military, data=samp))
model.14b <- summary(lm(policy ~ arm + age + gender + income + education + religious_id  + military, data=samp.threatonly))
model.15b <- summary(lm(policy ~ arm + age + gender + income + education + religious_id  + military, data=samp.unthreatonly))


#hypothesis 3: group####

#difference in means, classification

ate.19 <- mean(samp$action[samp$outgrp == 1],na.rm=T) - mean(samp$action[samp$ingrp == 1],na.rm=T) #vs. ingroup
ate.20 <- mean(samp$action[samp$outgrp == 1], na.rm=T) - mean(samp$action[samp$outgrp == 0], na.rm = T) #vs. all groups 


#standard errors of ates

se.ate.19 <- sqrt(var(samp$action[samp$outgrp == 1],na.rm=T)/
                    nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$action))) + 
                    var(samp$action[samp$ingrp == 1],na.rm=T)/
                    nrow(subset(samp, samp$ingrp == 1 & !is.na(samp$action))))

se.ate.20 <- sqrt(var(samp$action[samp$outgrp == 1],na.rm=T)/
                    nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$action))) + 
                    var(samp$action[samp$outgrp == 0],na.rm=T)/
                    nrow(subset(samp, samp$outgrp == 0 & !is.na(samp$action))))

#randomization inference

seri.ate.20 <- ri(samp, samp$action, 1000, 1017, 273)
seri.ate.20.ci <- c(seri.ate.20[round(.025*nrow(seri.ate.20)),1], seri.ate.20[round(.975*nrow(seri.ate.20)), 1])


#difference in means, repressive response

ate.21 <- mean(samp$policy[samp$outgrp == 1],na.rm=T) - mean(samp$policy[samp$ingrp == 1],na.rm=T) #vs. ingroup
ate.22 <- mean(samp$policy[samp$outgrp == 1], na.rm=T) - mean(samp$policy[samp$outgrp == 0], na.rm = T) #vs. all groups


#standard errors of ates

se.ate.21 <- sqrt(var(samp$policy[samp$outgrp == 1],na.rm=T)/
                    nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$policy))) + 
                    var(samp$policy[samp$ingrp == 1],na.rm=T)/
                    nrow(subset(samp, samp$ingrp == 1 & !is.na(samp$policy))))

se.ate.22 <- sqrt(var(samp$policy[samp$outgrp == 1],na.rm=T)/
                    nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$policy))) + 
                    var(samp$policy[samp$outgrp == 0],na.rm=T)/
                    nrow(subset(samp, samp$outgrp == 0 & !is.na(samp$policy))))

#randomization inference

seri.ate.22 <- ri(samp, samp$policy, 1000, 1017, 273)
seri.ate.22.ci <- c(seri.ate.22[round(.025*nrow(seri.ate.22)),1], seri.ate.22[round(.975*nrow(seri.ate.22)), 1])

#regression

samp.inoutonly <- subset(samp, samp$nogrp != 1)

model.19 <- summary(lm(action ~ outgrp, data=samp.inoutonly))
model.20 <- summary(lm(action ~ outgrp, data=samp))

model.21 <- summary(lm(policy ~ outgrp, data=samp.inoutonly))
model.22 <- summary(lm(policy ~ outgrp, data=samp))

#with covariates

model.19b <- summary(lm(action ~ outgrp + age + gender + income + education + religious_id  + military, data=samp.inoutonly))
model.20b <- summary(lm(action ~ outgrp + age + gender + income + education + religious_id  + military, data=samp))

model.21b <- summary(lm(policy ~ outgrp + age + gender + income + education + religious_id  + military, data=samp.inoutonly))
model.22b <- summary(lm(policy ~ outgrp + age + gender + income + education + religious_id  + military, data=samp))


#with interaction

model.19int <- summary(lm(action ~ outgrp*ideology, data=samp.inoutonly))
model.20int <- summary(lm(action ~ outgrp*ideology, data=samp))

model.21int <- summary(lm(policy ~ outgrp*ideology, data=samp.inoutonly))
model.22int <- summary(lm(policy ~ outgrp*ideology, data=samp))

#with covariates

model.19intb <- summary(lm(action ~ outgrp*ideology + age + gender + income + education + religious_id  + military, data=samp.inoutonly))
model.20intb <- summary(lm(action ~ outgrp*ideology+ age + gender + income + education + religious_id  + military, data=samp))

model.21intb <- summary(lm(policy ~ outgrp*ideology + age + gender + income + education + religious_id  + military, data=samp.inoutonly))
model.22intb <- summary(lm(policy ~ outgrp*ideology + age + gender + income + education + religious_id  + military, data=samp))


#hypothesis 3b: group and location####

#with interaction

model.25int <- summary(lm(action ~ outgrp*wbjeru, data=samp.inoutonly))
model.26int <- summary(lm(action ~ outgrp*wbjeru, data=samp))

model.27int <- summary(lm(policy ~ outgrp*wbjeru, data=samp.inoutonly))
model.28int <- summary(lm(policy ~ outgrp*wbjeru, data=samp))

#with covariates

model.25intb <- summary(lm(action ~ outgrp*wbjeru + age + gender + income + education + religious_id  + military + ideology, data=samp.inoutonly))
model.26intb <- summary(lm(action ~ outgrp*wbjeru+ age + gender + income + education + religious_id  + military + ideology, data=samp))

model.27intb <- summary(lm(policy ~ outgrp*wbjeru + age + gender + income + education + religious_id  + military + ideology, data=samp.inoutonly))
model.28intb <- summary(lm(policy ~ outgrp*wbjeru + age + gender + income + education + religious_id  + military + ideology, data=samp))


#hypothesis 3 extra: interaction of threat with ideology among subsamples####

samp.palestonly <- subset(samp, samp$outgrp == 1)
samp.israelonly <- subset(samp, samp$ingrp == 1)

#heterogeneous effects of threat by ideology
model.31int <- summary(lm(action ~ threat*ideology, data=samp.palestonly))
model.32int <- summary(lm(action ~ threat*ideology, data=samp.israelonly))

model.33int <- summary(lm(policy ~ threat*ideology, data=samp.palestonly))
model.34int <- summary(lm(policy ~ threat*ideology, data=samp.israelonly))


#heterogenous effects of arms by ideology
model.37int <- summary(lm(action ~ arm*ideology, data=samp.palestonly))
model.38int <- summary(lm(action ~ arm*ideology, data=samp.israelonly))

model.39int <- summary(lm(policy ~ arm*ideology, data=samp.palestonly))
model.40int <- summary(lm(policy ~ arm*ideology, data=samp.israelonly))

#4 plot results####

#hypothesis 1&2: threat and use of arms#### 
#plot coefficient estimates for change in perceptions based on threat

plot <- matrix(NA,6,3)

plot[,1] <- c(ate.1 - 1.96*se.ate.1, ate.2 - 1.96*se.ate.2, ate.3 - 1.96*se.ate.3,ate.10 - 1.96*se.ate.10, ate.11 - 1.96*se.ate.11, ate.12 - 1.96*se.ate.12)
plot[,2] <- c(ate.1,ate.2,ate.3,ate.10,ate.11,ate.12)
plot[,3] <- c(ate.1 + 1.96*se.ate.1, ate.2 + 1.96*se.ate.2, ate.3 + 1.96*se.ate.3,ate.10 + 1.96*se.ate.10, ate.11 + 1.96*se.ate.11, ate.12 + 1.96*se.ate.12)


plot<- as.data.frame(plot)
plot$names <- c("Threat","Threat (Armed Only)","Threat (Unarmed Only)", "Arms","Arms (Threat Only)","Arms (No Threat Only)")
plot$names <- as.factor(plot$names)
plot$names <- factor(plot$names, levels = plot$names[7 - row(as.matrix(plot$names))])
colnames(plot) <- c("lower","ptest","upper","names")

cfs1 <- ggplot(data = plot, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-0.25,0.75) + 
  xlab("Change in Proportion Violence") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50","gray50"))

plot2 <- matrix(NA,6,3)

plot2[,1] <- c(ate.4 - 1.96*se.ate.4, ate.5 - 1.96*se.ate.5, ate.6 - 1.96*se.ate.6, ate.13 - 1.96*se.ate.13, ate.14 - 1.96*se.ate.14, ate.15 - 1.96*ate.15)
plot2[,2] <- c(ate.4,ate.5,ate.6,ate.13,ate.14,ate.15)
plot2[,3] <- c(ate.4 + 1.96*se.ate.4, ate.5 + 1.96*se.ate.5, ate.6 + 1.96*se.ate.6, ate.13 + 1.96*se.ate.13, ate.14 + 1.96*se.ate.14, ate.15 + 1.96*ate.15)

plot2<- as.data.frame(plot2)
plot2$names <-  c("Threat","Threat (Armed Only)","Threat (Unarmed Only)", "Arms","Arms (Threat Only)","Arms (No Threat Only)")
plot2$names <- as.factor(plot2$names)
plot2$names <- factor(plot2$names, levels = plot2$names[7 - row(as.matrix(plot2$names))])
colnames(plot2) <- c("lower","ptest","upper","names")

cfs2 <- ggplot(data = plot2, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-0.25,0.75) + 
  xlab("Change in Repression Support") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50","gray50"))


#hypothesis 3: in-out group#### 

plot4 <- matrix(NA,2,3)

plot4[,1] <- c(ate.19 - 1.96*se.ate.19, ate.20 - 1.96*se.ate.20)
plot4[,2] <- c(ate.19, ate.20) 
plot4[,3] <- c(ate.19 + 1.96*se.ate.19, ate.20 + 1.96*se.ate.20)

plot4<- as.data.frame(plot4)
plot4$names <- c("Out vs. In","Out vs. Pooled")
plot4$names <- as.factor(plot4$names)
plot4$names <- factor(plot4$names, levels = plot4$names[3 - row(as.matrix(plot4$names))])
colnames(plot4) <- c("lower","ptest","upper","names")

cfs4 <- ggplot(data = plot4, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-0.5,0.75) + 
  xlab("Group Effects on Proportion Violence") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50"))

plot5 <- matrix(NA,2,3)

plot5[,1] <- c(ate.21 - 1.96*se.ate.21, ate.22 - 1.96*se.ate.22)
plot5[,2] <- c(ate.21, ate.22) 
plot5[,3] <- c(ate.21 + 1.96*se.ate.21, ate.22 + 1.96*se.ate.22)

plot5<- as.data.frame(plot5)
plot5$names <-  c("Out vs. In","Out vs. Pooled")
plot5$names <- as.factor(plot5$names)
plot5$names <- factor(plot5$names, levels = plot5$names[3 - row(as.matrix(plot5$names))])
colnames(plot5) <- c("lower","ptest","upper","names")

cfs5 <- ggplot(data = plot5, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-0.5,0.75) + 
  xlab("Group Effects on Repression Support") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50")) 


#interactions with right ideology####

plot7 <- matrix(NA,4,3)

plot7[,1] <- c(model.19int$coefficients[2,1] - 1.96*model.19int$coefficients[2,2], model.19int$coefficients[4,1] - 1.96*model.19int$coefficients[4,2], 
               model.20int$coefficients[2,1] - 1.96*model.20int$coefficients[2,2], model.20int$coefficients[4,1] - 1.96*model.20int$coefficients[4,2])
plot7[,2] <- c(model.19int$coefficients[2,1], model.19int$coefficients[4,1],model.20int$coefficients[2,1], model.20int$coefficients[4,1]) 
plot7[,3] <-  c(model.19int$coefficients[2,1] + 1.96*model.19int$coefficients[2,2], model.19int$coefficients[4,1] + 1.96*model.19int$coefficients[4,2], 
                model.20int$coefficients[2,1] + 1.96*model.20int$coefficients[2,2], model.20int$coefficients[4,1] + 1.96*model.20int$coefficients[4,2])

plot7<- as.data.frame(plot7)
plot7$names <- c("Palestinian vs. Settler","Pal vs. Set * Right","Palestinian vs. All", "Pal vs. All * Right")
plot7$names <- as.factor(plot7$names)
plot7$names <- factor(plot7$names, levels = plot7$names[5 - row(as.matrix(plot7$names))])
colnames(plot7) <- c("lower","ptest","upper","names")

cfs7 <- ggplot(data = plot7, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-1.25,0.75) + 
  xlab("Group Effects on Proportion Violence") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50"))

plot8 <- matrix(NA,4,3)

plot8[,1] <- c(model.21int$coefficients[2,1] - 1.96*model.21int$coefficients[2,2], model.21int$coefficients[4,1] - 1.96*model.21int$coefficients[4,2], 
               model.22int$coefficients[2,1] - 1.96*model.22int$coefficients[2,2], model.22int$coefficients[4,1] - 1.96*model.22int$coefficients[4,2])
plot8[,2] <- c(model.21int$coefficients[2,1], model.21int$coefficients[4,1],model.22int$coefficients[2,1], model.22int$coefficients[4,1]) 
plot8[,3] <-  c(model.21int$coefficients[2,1] + 1.96*model.21int$coefficients[2,2], model.21int$coefficients[4,1] + 1.96*model.21int$coefficients[4,2], 
                model.22int$coefficients[2,1] + 1.96*model.22int$coefficients[2,2], model.22int$coefficients[4,1] + 1.96*model.22int$coefficients[4,2])

plot8<- as.data.frame(plot8)
plot8$names <- c("Palestinian vs. Settler","Pal vs. Set * Right","Palestinian vs. All", "Pal vs. All * Right")
plot8$names <- as.factor(plot8$names)
plot8$names <- factor(plot8$names, levels = plot8$names[5 - row(as.matrix(plot8$names))])
colnames(plot8) <- c("lower","ptest","upper","names")

cfs8 <- ggplot(data = plot8, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-1.25,0.75) + 
  xlab("Group Effects on Repression Support") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50")) 



#interactions with location####

plot10 <- matrix(NA,4,3)

plot10[,1] <- c(model.25int$coefficients[2,1] - 1.96*model.25int$coefficients[2,2], model.25int$coefficients[4,1] - 1.96*model.25int$coefficients[4,2], 
                model.26int$coefficients[2,1] - 1.96*model.26int$coefficients[2,2], model.26int$coefficients[4,1] - 1.96*model.26int$coefficients[4,2])
plot10[,2] <- c(model.25int$coefficients[2,1], model.25int$coefficients[4,1],model.26int$coefficients[2,1], model.26int$coefficients[4,1]) 
plot10[,3] <-  c(model.25int$coefficients[2,1] + 1.96*model.25int$coefficients[2,2], model.25int$coefficients[4,1] + 1.96*model.25int$coefficients[4,2], 
                 model.26int$coefficients[2,1] + 1.96*model.26int$coefficients[2,2], model.26int$coefficients[4,1] + 1.96*model.26int$coefficients[4,2])

plot10<- as.data.frame(plot10)
plot10$names <- c("Palestinian vs. Settler","Pal vs. Set * WB/Jerusalem","Palestinian vs. All", "Pal vs. All * WB/Jerusalem")
plot10$names <- as.factor(plot10$names)
plot10$names <- factor(plot10$names, levels = plot10$names[5 - row(as.matrix(plot10$names))])
colnames(plot10) <- c("lower","ptest","upper","names")

cfs10 <- ggplot(data = plot10, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-0.5,0.75) + 
  xlab("Group Effects on Proportion Violence") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50"))

plot11 <- matrix(NA,4,3)

plot11[,1] <- c(model.27int$coefficients[2,1] - 1.96*model.27int$coefficients[2,2], model.27int$coefficients[4,1] - 1.96*model.27int$coefficients[4,2], 
                model.28int$coefficients[2,1] - 1.96*model.28int$coefficients[2,2], model.28int$coefficients[4,1] - 1.96*model.28int$coefficients[4,2])
plot11[,2] <- c(model.27int$coefficients[2,1], model.27int$coefficients[4,1],model.28int$coefficients[2,1], model.28int$coefficients[4,1]) 
plot11[,3] <-  c(model.27int$coefficients[2,1] + 1.96*model.27int$coefficients[2,2], model.27int$coefficients[4,1] + 1.96*model.27int$coefficients[4,2], 
                 model.28int$coefficients[2,1] + 1.96*model.28int$coefficients[2,2], model.28int$coefficients[4,1] + 1.96*model.28int$coefficients[4,2])

plot11<- as.data.frame(plot11)
plot11$names <- c("Palestinian vs. Settler","Pal vs. Set * WB/Jerusalem","Palestinian vs. All", "Pal vs. All * WB/Jerusalem")
plot11$names <- as.factor(plot11$names)
plot11$names <- factor(plot11$names, levels = plot11$names[5 - row(as.matrix(plot11$names))])
colnames(plot11) <- c("lower","ptest","upper","names")

cfs11 <- ggplot(data = plot11, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-0.5,0.75) + 
  xlab("Group Effects on Repression Support") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50")) 



library(gridExtra)


#############################################
################Figure 1 and 2 in Main Text
############################################
grid.arrange(cfs1,cfs2,ncol=2)
grid.arrange(cfs4,cfs5,ncol=2)

# plots correspond to Tables 8 and 9 in Appendix
# grid.arrange(cfs7,cfs8,ncol=2)
# grid.arrange(cfs10,cfs11,ncol=2)

#5 tables####

library(stargazer)

#reformat for stargazer

#############################
## Table 4 in Appendix
#############################

model.1 <- lm(action ~ threat, data=samp)
model.2 <- lm(action ~ threat, data=samp.armedonly)
model.3 <- lm(action ~ threat, data=samp.unarmedonly)

model.1b <- lm(action ~ threat + age + gender + income + education + religious_id  + military, data=samp)
model.2b <- lm(action ~ threat + age + gender + income + education + religious_id  + military, data=samp.armedonly)
model.3b <- lm(action ~ threat + age + gender + income + education + religious_id  + military, data=samp.unarmedonly)

stargazer(model.1, model.2, model.3, model.1b, model.2b, model.3b,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Perception of Violence"),
          covariate.labels=c("Threat of Harm","Age","Gender (Male)", "Gender (Other)","Income","Education","Religious","Military"), 
          out="exp_models_israel_1.txt")


#############################
## Table 5 in Appendix
#############################


model.4 <- lm(policy ~ threat, data=samp)
model.5 <- lm(policy ~ threat, data=samp.armedonly)
model.6 <- lm(policy ~ threat, data=samp.unarmedonly)

model.4b <- lm(policy ~ threat + age + gender + income + education + religious_id  + military, data=samp)
model.5b <- lm(policy ~ threat + age + gender + income + education + religious_id  + military, data=samp.armedonly)
model.6b <- lm(policy ~ threat + age + gender + income + education + religious_id  + military, data=samp.unarmedonly)

stargazer(model.4, model.5, model.6, model.4b, model.5b, model.6b,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Support for Repression"),
          covariate.labels=c("Threat of Harm","Age","Gender (Male)", "Gender (Other)","Income","Education","Religious","Military"), 
          out="exp_models_israel_2.txt")


#############################
## Table 6 in Appendix
#############################


model.10 <- lm(action ~ arm, data=samp)
model.11 <- lm(action ~ arm, data=samp.threatonly)
model.12 <- lm(action ~ arm, data=samp.unthreatonly)

model.10b <- lm(action ~ arm + age + gender + income + education + religious_id  + military, data=samp)
model.11b <- lm(action ~ arm + age + gender + income + education + religious_id  + military, data=samp.threatonly)
model.12b <- lm(action ~ arm + age + gender + income + education + religious_id  + military, data=samp.unthreatonly)

stargazer(model.10, model.11, model.12, model.10b, model.11b, model.12b,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Perception of Violence"),
          covariate.labels=c("Use of Arms","Age","Gender (Male)", "Gender (Other)","Income","Education","Religious","Military"), 
          out="exp_models_israel_4.txt")


#############################
## Table 7 in Appendix
#############################

model.13 <- lm(policy ~ arm, data=samp)
model.14 <- lm(policy ~ arm, data=samp.threatonly)
model.15 <- lm(policy ~ arm, data=samp.unthreatonly)

model.13b <- lm(policy ~ arm + age + gender + income + education + religious_id  + military, data=samp)
model.14b <- lm(policy ~ arm + age + gender + income + education + religious_id  + military, data=samp.threatonly)
model.15b <- lm(policy ~ arm + age + gender + income + education + religious_id  + military, data=samp.unthreatonly)

stargazer(model.13, model.14, model.15, model.13b, model.14b, model.15b,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Support for Repression"),
          covariate.labels=c("Use of Arms","Age","Gender (Male)", "Gender (Other)","Income","Education","Religious","Military"), 
          out="exp_models_israel_5.txt")


#############################
## Table 8 in Appendix
#############################


model.19int <- lm(action ~ outgrp*ideology, data=samp.inoutonly)
model.20int <- lm(action ~ outgrp*ideology, data=samp)
model.19intb <- lm(action ~ outgrp*ideology + age + gender + income + education + religious_id  + military, data=samp.inoutonly)
model.20intb <- lm(action ~ outgrp*ideology + age + gender + income + education + religious_id  + military, data=samp)

stargazer(model.19int, model.20int, model.19intb, model.20intb,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Perception of Violence"),
          covariate.labels=c("Palestinian","Right Ideology","Age","Gender (Male)", "Gender (Other)","Income","Education","Religious","Military","Palest*Right"), 
          out="exp_models_israel_7.txt")


#############################
## Table 9 in Appendix
#############################

model.21int <- lm(policy ~ outgrp*ideology, data=samp.inoutonly)
model.22int <- lm(policy ~ outgrp*ideology, data=samp)
model.21intb <- lm(policy ~ outgrp*ideology + age + gender + income + education + religious_id  + military, data=samp.inoutonly)
model.22intb <- lm(policy ~ outgrp*ideology+ age + gender + income + education + religious_id  + military, data=samp)

stargazer(model.21int, model.22int, model.21intb, model.22intb,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Support for Repression"),
          covariate.labels=c("Palestinian","Right Ideology","Age","Gender (Male)","Gender (Other)","Income","Education","Religious","Military","Palest*Right"), 
          out="exp_models_israel_8.txt")


#############################
## Table 10 in Appendix
#############################


model.43 <- lm(crimehatemotive ~ threat, data=samp)
model.44 <- lm(crimehatemotive ~ threat, data=samp.armedonly)
model.45 <- lm(crimehatemotive ~ threat, data=samp.unarmedonly)

model.43b <- lm(crimehatemotive ~ threat + age + gender + income + education + religious_id  + military, data=samp)
model.44b <- lm(crimehatemotive ~ threat + age + gender + income + education + religious_id  + military, data=samp.armedonly)
model.45b <- lm(crimehatemotive ~ threat + age + gender + income + education + religious_id  + military, data=samp.unarmedonly)

stargazer(model.43, model.44, model.45, model.43b, model.44b, model.45b,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Assign as Crime or Hate"),
          covariate.labels=c("Threat of Harm","Age","Gender (Male)", "Gender (Other)","Income","Education","Religious","Military"), 
          out="exp_models_israel_13.txt")


#############################
## Table 11 in Appendix
#############################


model.46 <- lm(crimehatemotive ~ arm, data=samp)
model.47 <- lm(crimehatemotive ~ arm, data=samp.threatonly)
model.48 <- lm(crimehatemotive ~ arm, data=samp.unthreatonly)

model.46b <- lm(crimehatemotive ~ arm + age + gender + income + education + religious_id  + military, data=samp)
model.47b <- lm(crimehatemotive ~ arm + age + gender + income + education + religious_id  + military, data=samp.threatonly)
model.48b <- lm(crimehatemotive ~ arm + age + gender + income + education + religious_id  + military, data=samp.unthreatonly)

stargazer(model.46, model.47, model.48, model.46b, model.47b, model.48b,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Assign as Crime or Hate"),
          covariate.labels=c("Use of Arms","Age","Gender (Male)", "Gender (Other)","Income","Education","Religious","Military"), 
          out="exp_models_israel_14.txt")


#############################
## Table 12 in Appendix
#############################

model.49 <- lm(crimehatemotive ~ outgrp, data=samp)
model.50 <- lm(crimehatemotive ~ outgrp, data=samp.inoutonly)

model.49b <- lm(crimehatemotive ~ outgrp + age + gender + income + education + religious_id  + military, data=samp)
model.50b <- lm(crimehatemotive ~ outgrp + age + gender + income + education + religious_id  + military, data=samp.inoutonly)

stargazer(model.49, model.50, model.49b, model.50b,
          type="latex", 
          column.separate=NULL,
          dep.var.labels=c("Assign as Crime or Hate"),
          covariate.labels=c("Palestinian","Age","Gender (Male)","Gender (Other)", "Income","Education","Religious","Military"), 
          out="exp_models_israel_15.txt")


################################################
#balance tests####
# Figure 3 in Appendix
####################################################

names <- c("threat","arm","outgrp","education","income","ideology","age","religious_id","military")
samp.bal <- samp[names]
samp.bal <- na.omit(samp.bal)

#threat treatment

library(ggplot2)
library(ggthemes)

bal <- matrix(NA,6,3)

bal[,2] <- c((mean(samp.bal$education[samp.bal$threat==1],na.rm=T)-mean(samp.bal$education[samp.bal$threat==0],na.rm=T))/sd(samp.bal$education[samp.bal$threat==1],na.rm=T),
             (mean(samp.bal$income[samp.bal$threat==1],na.rm=T)-mean(samp.bal$income[samp.bal$threat==0],na.rm=T))/sd(samp.bal$income[samp.bal$threat==1],na.rm=T),
             (mean(samp.bal$religious_id[samp.bal$threat==1],na.rm=T)-mean(samp.bal$religious_id[samp.bal$threat==0],na.rm=T))/sd(samp.bal$religious_id[samp.bal$threat==1],na.rm=T),         
             (mean(samp.bal$military[samp.bal$threat==1],na.rm=T)-mean(samp.bal$military[samp.bal$threat==0],na.rm=T))/sd(samp.bal$military[samp.bal$threat==1],na.rm=T),
             (mean(samp.bal$ideology[samp.bal$threat==1],na.rm=T)-mean(samp.bal$ideology[samp.bal$threat==0],na.rm=T))/sd(samp.bal$ideology[samp.bal$threat==1],na.rm=T),
             (mean(samp.bal$age[samp.bal$threat==1],na.rm=T)-mean(samp.bal$age[samp.bal$threat==0],na.rm=T))/sd(samp.bal$age[samp.bal$threat==1],na.rm=T))

treatnum <- sum(samp.bal$threat==1)
controlnum <- sum(samp.bal$threat==0)

bal[,1] <- c(bal[1,2] - 1.96*sqrt(var(samp.bal$education[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$threat==1],na.rm=T),
             bal[2,2] - 1.96*sqrt(var(samp.bal$income[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$threat==1],na.rm=T),
             bal[3,2] - 1.96*sqrt(var(samp.bal$religious_id[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$religious_id[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$religious_id[samp.bal$threat==1],na.rm=T),
             bal[4,2] - 1.96*sqrt(var(samp.bal$military[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$military[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$military[samp.bal$threat==1],na.rm=T),
             bal[5,2] - 1.96*sqrt(var(samp.bal$ideology[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$threat==1],na.rm=T),
             bal[6,2] - 1.96*sqrt(var(samp.bal$age[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$threat==1],na.rm=T))


bal[,3] <- c(bal[1,2] + 1.96*sqrt(var(samp.bal$education[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$threat==1],na.rm=T),
             bal[2,2] + 1.96*sqrt(var(samp.bal$income[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$threat==1],na.rm=T),
             bal[3,2] + 1.96*sqrt(var(samp.bal$religious_id[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$religious_id[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$religious_id[samp.bal$threat==1],na.rm=T),
             bal[4,2] + 1.96*sqrt(var(samp.bal$military[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$military[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$military[samp.bal$threat==1],na.rm=T),
             bal[5,2] + 1.96*sqrt(var(samp.bal$ideology[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$threat==1],na.rm=T),
             bal[6,2] + 1.96*sqrt(var(samp.bal$age[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$threat==1],na.rm=T))

bal <- as.data.frame(bal)
bal$names <- c("Education","Income","Religious ID","Military Service","Ideology","Age")
colnames(bal) <- c("lower","diff","upper","names")

balplot <- ggplot(data = bal, aes(x = diff, y = names, color = names)) +
  geom_point(size=1.25) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1) +
  xlim(-0.4,0.4) + 
  xlab("Threat Treatments") +
  ylab("Covariates") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50","gray50"))

#arms treatment

bal2 <- matrix(NA,6,3)

bal2[,2] <- c((mean(samp.bal$education[samp.bal$arm==1],na.rm=T)-mean(samp.bal$education[samp.bal$arm==0],na.rm=T))/sd(samp.bal$education[samp.bal$arm==1],na.rm=T),
              (mean(samp.bal$income[samp.bal$arm==1],na.rm=T)-mean(samp.bal$income[samp.bal$arm==0],na.rm=T))/sd(samp.bal$income[samp.bal$arm==1],na.rm=T),
              (mean(samp.bal$religious_id[samp.bal$arm==1],na.rm=T)-mean(samp.bal$religious_id[samp.bal$arm==0],na.rm=T))/sd(samp.bal$religious_id[samp.bal$arm==1],na.rm=T),      
              (mean(samp.bal$religious_id[samp.bal$arm==1],na.rm=T)-mean(samp.bal$religious_id[samp.bal$arm==0],na.rm=T))/sd(samp.bal$religious_id[samp.bal$arm==1],na.rm=T),         
              (mean(samp.bal$ideology[samp.bal$arm==1],na.rm=T)-mean(samp.bal$ideology[samp.bal$arm==0],na.rm=T))/sd(samp.bal$ideology[samp.bal$arm==1],na.rm=T),
              (mean(samp.bal$age[samp.bal$arm==1],na.rm=T)-mean(samp.bal$age[samp.bal$arm==0],na.rm=T))/sd(samp.bal$age[samp.bal$arm==1],na.rm=T))

treatnum <- sum(samp.bal$arm==1)
controlnum <- sum(samp.bal$arm==0)

bal2[,1] <- c(bal2[1,2] - 1.96*sqrt(var(samp.bal$education[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$arm==1],na.rm=T),
              bal2[2,2] - 1.96*sqrt(var(samp.bal$income[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$arm==1],na.rm=T),
              bal2[3,2] - 1.96*sqrt(var(samp.bal$religious_id[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$religious_id[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$religious_id[samp.bal$arm==1],na.rm=T),
              bal2[4,2] - 1.96*sqrt(var(samp.bal$military[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$military[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$military[samp.bal$arm==1],na.rm=T),
              bal2[5,2] - 1.96*sqrt(var(samp.bal$ideology[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$arm==1],na.rm=T),
              bal2[6,2] - 1.96*sqrt(var(samp.bal$age[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$arm==1],na.rm=T))


bal2[,3] <- c(bal2[1,2] + 1.96*sqrt(var(samp.bal$education[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$arm==1],na.rm=T),
              bal2[2,2] + 1.96*sqrt(var(samp.bal$income[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$arm==1],na.rm=T),
              bal2[3,2] + 1.96*sqrt(var(samp.bal$religious_id[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$religious_id[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$religious_id[samp.bal$arm==1],na.rm=T),
              bal2[4,2] + 1.96*sqrt(var(samp.bal$military[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$military[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$military[samp.bal$arm==1],na.rm=T),
              bal2[5,2] + 1.96*sqrt(var(samp.bal$ideology[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$arm==1],na.rm=T),
              bal2[6,2] + 1.96*sqrt(var(samp.bal$age[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$arm==1],na.rm=T))

bal2 <- as.data.frame(bal2)
bal2$names <- c("Education","Income","Religious ID","Military Service","Ideology","Age")
colnames(bal2) <- c("lower","diff","upper","names")

balplot2 <- ggplot(data = bal2, aes(x = diff, y = names, color = names)) +
  geom_point(size=1.25) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1) +
  xlim(-0.4,0.4) + 
  xlab("Arms Treatment") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none")+
  scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50","gray50"))


#outgroup treatment

bal3 <- matrix(NA,6,3)

bal3[,2] <- c((mean(samp.bal$education[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$education[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$education[samp.bal$outgrp==1],na.rm=T),
              (mean(samp.bal$income[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$income[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$income[samp.bal$outgrp==1],na.rm=T),
              (mean(samp.bal$religious_id[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$religious_id[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$religious_id[samp.bal$outgrp==1],na.rm=T), 
              (mean(samp.bal$military[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$military[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$military[samp.bal$outgrp==1],na.rm=T),         
              (mean(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$ideology[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T),
              (mean(samp.bal$age[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$age[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$age[samp.bal$outgrp==1],na.rm=T))

treatnum <- sum(samp.bal$outgrp==1)
controlnum <- sum(samp.bal$outgrp==0)

bal3[,1] <- c(bal3[1,2] - 1.96*sqrt(var(samp.bal$education[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$outgrp==1],na.rm=T),
              bal3[2,2] - 1.96*sqrt(var(samp.bal$income[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$outgrp==1],na.rm=T),
              bal3[3,2] - 1.96*sqrt(var(samp.bal$religious_id[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$religious_id[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$religious_id[samp.bal$outgrp==1],na.rm=T),
              bal3[4,2] - 1.96*sqrt(var(samp.bal$military[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$military[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$religious_id[samp.bal$outgrp==1],na.rm=T),
              bal3[5,2] - 1.96*sqrt(var(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T),
              bal3[6,2] - 1.96*sqrt(var(samp.bal$age[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$outgrp==1],na.rm=T))


bal3[,3] <- c(bal3[1,2] + 1.96*sqrt(var(samp.bal$education[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$outgrp==1],na.rm=T),
              bal3[2,2] + 1.96*sqrt(var(samp.bal$income[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$outgrp==1],na.rm=T),
              bal3[3,2] + 1.96*sqrt(var(samp.bal$religious_id[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$religious_id[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$religious_id[samp.bal$outgrp==1],na.rm=T),
              bal3[4,2] + 1.96*sqrt(var(samp.bal$military[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$military[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$military[samp.bal$outgrp==1],na.rm=T),
              bal3[5,2] + 1.96*sqrt(var(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T),
              bal3[6,2] + 1.96*sqrt(var(samp.bal$age[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$outgrp==1],na.rm=T))

bal3 <- as.data.frame(bal3)
bal3$names <- c("Education","Income","Religious ID","Military Service","Ideology","Age")
colnames(bal3) <- c("lower","diff","upper","names")

balplot3 <- ggplot(data = bal3, aes(x = diff, y = names, color = names)) +
  geom_point(size=1.25) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1) +
  xlim(-0.41,0.4) + 
  xlab("Out-Group Treatment") +
  ylab(" ") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50","gray50"))

library(gridExtra)

################################################
#balance tests####
# Figure 3 in Appendix
####################################################

grid.arrange(balplot,balplot2,balplot3, ncol=3)

#attrition####

######################
## attrition bias check
#################

#MIPO?
attrition <- summary(lm(manipulation ~ factor(treatment) + gender + age + education + income + party + religious_id, data=samp))

#test with randomizataion inference according to Gerber and Green
ri.lm.multiple <- function(df, outcome, covariate1, covariate2, covariate3, covariate4, covariate5, covariate6, reps, n, m, num){
  
  library(randomizr)
  #create objects
  Fs <- vector(mode="numeric")
  assns <- as.data.frame(matrix(NA,nrow=nrow(df),ncol=reps))
  dfs <- list()
  
  #generate random assignments
  for(i in 1:reps){
    assns[,i] <- complete_ra(N = n, prob_each = m, num_arms = num)}
  
  #remove duplicate assignments
  n=reps
  assns <- assns[!duplicated(lapply(assns, function(k) head(k, n)))]
  
  #take F statistic from regression and store
  for(i in 1:ncol(assns)){dfs[[i]] <- cbind(outcome, assns[,i], covariate1,covariate2, covariate3, covariate4, covariate5, covariate6) 
  Fs[i] <- summary(lm(dfs[[i]][,1] ~ factor(dfs[[i]][,2]) + dfs[[i]][,3]+ dfs[[i]][,4]+ dfs[[i]][,5]+ 
                        dfs[[i]][,6]+ dfs[[i]][,7] + dfs[[i]][,8]))$fstatistic[1]}
  Fs <- as.data.frame(sort(Fs))
  
  #output
  return(Fs)
}

fstats <- ri.lm.multiple(samp, samp$manipulation, samp$gender, samp$age, samp$education, samp$income, samp$party, samp$religious_id,
                         1000,1017,rep(1/12,12),12)

sig <- sum(ifelse(fstats >= attrition$fstatistic[1],1,0))


#####################################################################################
############END OF ISRAEL SURVEY
############# MAKE SURE TO CLEAR VARIABLE NAMES BETWEEN US AND ISRAEL SURVEY
############# CODE WILL BREAK OTHERWISE, DUE TO NAME REPITITIONS
#################################################################################


#############################
#US Survey
##############################
rm(list = ls())
cat("\014")
set.seed(30030)

samp <- read.csv("US_Surveydata_final.csv")


############################################
#### Main Hypothesis tests
#### Figures 3 and 4 in Main Text 
###########################################


#hypothesis 1: threat####

#difference in means, classification

ate.1 <- mean(samp$action[samp$threat == 1],na.rm=T) - mean(samp$action[samp$threat == 0],na.rm=T) #all
ate.2 <- mean(samp$action[samp$arm_conf == 1], na.rm=T) - mean(samp$action[samp$arm_enc == 1], na.rm = T) #within armed
ate.3 <- mean(samp$action[samp$un_conf == 1], na.rm=T) - mean(samp$action[samp$pr_nv == 1], na.rm = T) #within unarmed

  
  #standard errors of ates
  
  se.ate.1 <- sqrt(var(samp$action[samp$threat == 1],na.rm=T)/
                     nrow(subset(samp, samp$threat == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$threat == 0],na.rm=T)/
                     nrow(subset(samp, samp$threat == 0 & !is.na(samp$action))))
  
  se.ate.2 <- sqrt(var(samp$action[samp$arm_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$arm_enc == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$action))))
  
  se.ate.3 <- sqrt(var(samp$action[samp$un_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$pr_nv == 1],na.rm=T)/
                     nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$action))))
  
  #randomization inference
library(randomizr)  
  
  ri <- function(df, outcome, reps, n, m){
    
    library(randomizr)
    #create objects
    ATEs <- vector(mode="numeric")
    assns <- as.data.frame(matrix(NA,nrow=nrow(df),ncol=reps))
    dfs <- list()
    
    #generate random assignments
    for(i in 1:reps){
      assns[,i] <- complete_ra(N = n, m = m)}
    
    #remove duplicate assignments
    n=reps
    assns <- assns[!duplicated(lapply(assns, function(k) head(k, n)))]
    
    #calculate ATE for each assignment and store
    for(i in 1:ncol(assns)){dfs[[i]] <- cbind(outcome, assns[,i]) 
    ATEs[i] <- mean(dfs[[i]][,1][dfs[[i]][,2] == 1],na.rm=T) - mean(dfs[[i]][,1][dfs[[i]][,2] == 0],na.rm=T)}
    ATEs <- as.data.frame(sort(ATEs))
    
    #output
    return(ATEs)
  }
  
  seri.ate.1 <- ri(samp, samp$action, 1000, 1137, 530)
  seri.ate.1.ci <- c(seri.ate.1[round(.025*nrow(seri.ate.1)),1], seri.ate.1[round(.975*nrow(seri.ate.1)), 1])
  
#difference in means, repressive response
  
  ate.4 <- mean(samp$policy[samp$threat == 1],na.rm=T) - mean(samp$policy[samp$threat == 0],na.rm=T) #all
  ate.5 <- mean(samp$policy[samp$arm_conf == 1], na.rm=T) - mean(samp$policy[samp$arm_enc == 1], na.rm = T) #within armed
  ate.6 <- mean(samp$policy[samp$un_conf == 1], na.rm=T) - mean(samp$policy[samp$pr_nv == 1], na.rm = T) #within unarmed
  
  
  #standard errors of ates
  
  se.ate.4 <- sqrt(var(samp$policy[samp$threat == 1],na.rm=T)/
                     nrow(subset(samp, samp$threat == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$threat == 0],na.rm=T)/
                     nrow(subset(samp, samp$threat == 0 & !is.na(samp$policy))))
  
  se.ate.5 <- sqrt(var(samp$policy[samp$arm_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$arm_enc == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$policy))))
  
  se.ate.6 <- sqrt(var(samp$policy[samp$un_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$pr_nv == 1],na.rm=T)/
                     nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$policy))))
  
  #randomization inference
  
  seri.ate.4 <- ri(samp, samp$policy, 1000, 1137, 530)
  seri.ate.4.ci <- c(seri.ate.4[round(.025*nrow(seri.ate.4)),1], seri.ate.4[round(.975*nrow(seri.ate.4)), 1])

#regression
    
  samp.armedonly <- subset(samp, samp$arm == 1)
  samp.unarmedonly <- subset(samp, samp$arm == 0)

  model.1 <- summary(lm(action ~ threat, data=samp))
  model.2 <- summary(lm(action ~ threat, data=samp.armedonly))
  model.3 <- summary(lm(action ~ threat, data=samp.unarmedonly))
  
  model.4 <- summary(lm(policy ~ threat, data=samp))
  model.5 <- summary(lm(policy ~ threat, data=samp.armedonly))
  model.6 <- summary(lm(policy ~ threat, data=samp.unarmedonly))
  
  #with covariates
  
  model.1b <- summary(lm(action ~ threat + age + gender + income + education + pid7, data=samp))
  model.2b <- summary(lm(action ~ threat + age + gender + income + education + pid7, data=samp.armedonly))
  model.3b <- summary(lm(action ~ threat + age + gender + income + education + pid7, data=samp.unarmedonly))
  
  model.4b <- summary(lm(policy ~ threat + age + gender + income + education + pid7, data=samp))
  model.5b <- summary(lm(policy ~ threat + age + gender + income + education + pid7, data=samp.armedonly))
  model.6b <- summary(lm(policy ~ threat + age + gender + income + education + pid7, data=samp.unarmedonly))
  

#hypothesis 2: arms####

  #difference in means, classification
  
  ate.10 <- mean(samp$action[samp$arm == 1],na.rm=T) - mean(samp$action[samp$arm == 0],na.rm=T) #all
  ate.11 <- mean(samp$action[samp$arm_conf == 1], na.rm=T) - mean(samp$action[samp$un_conf == 1], na.rm = T) #within threatening
  ate.12 <- mean(samp$action[samp$arm_enc == 1], na.rm=T) - mean(samp$action[samp$pr_nv == 1], na.rm = T) #within unthreatening
  
  #standard errors of ates
  
  se.ate.10 <- sqrt(var(samp$action[samp$arm == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$arm == 0],na.rm=T)/
                     nrow(subset(samp, samp$arm == 0 & !is.na(samp$action))))
  
  se.ate.11 <- sqrt(var(samp$action[samp$arm_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$un_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$action))))
  
  se.ate.12 <- sqrt(var(samp$action[samp$arm_enc == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$pr_nv == 1],na.rm=T)/
                     nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$action))))
  
  #randomization inference
  
  seri.ate.10 <- ri(samp, samp$action, 1000, 1137, 533)
  seri.ate.10.ci <- c(seri.ate.10[round(.025*nrow(seri.ate.10)),1], seri.ate.10[round(.975*nrow(seri.ate.10)), 1])
  
#difference in means, repressive response
  
  ate.13 <- mean(samp$policy[samp$arm == 1],na.rm=T) - mean(samp$policy[samp$arm == 0],na.rm=T) #all
  ate.14 <- mean(samp$policy[samp$arm_conf == 1], na.rm=T) - mean(samp$policy[samp$un_conf == 1], na.rm = T) #within threatening
  ate.15 <- mean(samp$policy[samp$arm_enc == 1], na.rm=T) - mean(samp$policy[samp$pr_nv == 1], na.rm = T) #within unthreatening
  
  
  #standard errors of ates
  
  se.ate.13 <- sqrt(var(samp$policy[samp$arm == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$arm == 0],na.rm=T)/
                     nrow(subset(samp, samp$arm == 0 & !is.na(samp$policy))))
  
  se.ate.14 <- sqrt(var(samp$policy[samp$arm_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_conf == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$un_conf == 1],na.rm=T)/
                     nrow(subset(samp, samp$un_conf == 1 & !is.na(samp$policy))))
  
  se.ate.15 <- sqrt(var(samp$policy[samp$arm_enc == 1],na.rm=T)/
                     nrow(subset(samp, samp$arm_enc == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$pr_nv == 1],na.rm=T)/
                     nrow(subset(samp, samp$pr_nv == 1 & !is.na(samp$policy))))
  
  #randomization inference
  
  seri.ate.13 <- ri(samp, samp$policy, 1000, 1137, 533)
  seri.ate.13.ci <- c(seri.ate.13[round(.025*nrow(seri.ate.13)),1], seri.ate.13[round(.975*nrow(seri.ate.13)), 1])
  
  
#regression
  
  samp.threatonly <- subset(samp, samp$threat == 1)
  samp.unthreatonly <- subset(samp, samp$threat == 0)
  
  model.10 <- summary(lm(action ~ arm, data=samp))
  model.11 <- summary(lm(action ~ arm, data=samp.threatonly))
  model.12 <- summary(lm(action ~ arm, data=samp.unthreatonly))
  
  model.13 <- summary(lm(policy ~ arm, data=samp))
  model.14 <- summary(lm(policy ~ arm, data=samp.threatonly))
  model.15 <- summary(lm(policy ~ arm, data=samp.unthreatonly))
  
  model.10b <- summary(lm(action ~ arm + age + gender + income + education + pid7, data=samp))
  model.11b <- summary(lm(action ~ arm + age + gender + income + education + pid7, data=samp.threatonly))
  model.12b <- summary(lm(action ~ arm + age + gender + income + education + pid7, data=samp.unthreatonly))
  
  model.13b <- summary(lm(policy ~ arm + age + gender + income + education + pid7, data=samp))
  model.14b <- summary(lm(policy ~ arm + age + gender + income + education + pid7, data=samp.threatonly))
  model.15b <- summary(lm(policy ~ arm + age + gender + income + education + pid7, data=samp.unthreatonly))
  
#hypothesis 3: group####
  
  #difference in means, classification
  
  ate.19 <- mean(samp$action[samp$outgrp == 1],na.rm=T) - mean(samp$action[samp$ingrp == 1],na.rm=T) #vs. ingroup
  ate.20 <- mean(samp$action[samp$outgrp == 1], na.rm=T) - mean(samp$action[samp$outgrp == 0], na.rm = T) #vs. all groups 
  
  
  #standard errors of ates
  
  se.ate.19 <- sqrt(var(samp$action[samp$outgrp == 1],na.rm=T)/
                     nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$ingrp == 1],na.rm=T)/
                     nrow(subset(samp, samp$ingrp == 1 & !is.na(samp$action))))
  
  se.ate.20 <- sqrt(var(samp$action[samp$outgrp == 1],na.rm=T)/
                     nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$action))) + 
                     var(samp$action[samp$outgrp == 0],na.rm=T)/
                     nrow(subset(samp, samp$outgrp == 0 & !is.na(samp$action))))
  
  
  #difference in means, repressive response
  
  ate.21 <- mean(samp$policy[samp$outgrp == 1],na.rm=T) - mean(samp$policy[samp$ingrp == 1],na.rm=T) #vs. ingroup
  ate.22 <- mean(samp$policy[samp$outgrp == 1], na.rm=T) - mean(samp$policy[samp$outgrp == 0], na.rm = T) #vs. all groups
  
  
  #standard errors of ates
  
  se.ate.21 <- sqrt(var(samp$policy[samp$outgrp == 1],na.rm=T)/
                     nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$ingrp == 1],na.rm=T)/
                     nrow(subset(samp, samp$ingrp == 1 & !is.na(samp$policy))))
  
  se.ate.22 <- sqrt(var(samp$policy[samp$outgrp == 1],na.rm=T)/
                     nrow(subset(samp, samp$outgrp == 1 & !is.na(samp$policy))) + 
                     var(samp$policy[samp$outgrp == 0],na.rm=T)/
                     nrow(subset(samp, samp$outgrp == 0 & !is.na(samp$policy))))

  

  samp.inoutonly <- subset(samp, samp$nogrp != 1)
  
  model.19 <- summary(lm(action ~ outgrp, data=samp.inoutonly))
  model.20 <- summary(lm(action ~ outgrp, data=samp))

  model.21 <- summary(lm(policy ~ outgrp, data=samp.inoutonly))
  model.22 <- summary(lm(policy ~ outgrp, data=samp))
  
  #with covariates
  
  model.19b <- summary(lm(action ~ outgrp + age + gender + income + education + pid7, data=samp.inoutonly))
  model.20b <- summary(lm(action ~ outgrp + age + gender + income + education + pid7, data=samp))
  
  model.21b <- summary(lm(policy ~ outgrp + age + gender + income + education + pid7, data=samp.inoutonly))
  model.22b <- summary(lm(policy ~ outgrp + age + gender + income + education + pid7, data=samp))
  
  #with interaction
  
  model.19int <- summary(lm(action ~ outgrp*gop_id, data=samp.inoutonly))
  model.20int <- summary(lm(action ~ outgrp*gop_id, data=samp))
  
  model.21int <- summary(lm(policy ~ outgrp*gop_id, data=samp.inoutonly))
  model.22int <- summary(lm(policy ~ outgrp*gop_id, data=samp))
  
  
  #with covariates
  
  model.19intb <- summary(lm(action ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp.inoutonly))
  model.20intb <- summary(lm(action ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp))
  
  model.21intb <- summary(lm(policy ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp.inoutonly))
  model.22intb <- summary(lm(policy ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp))
  
  
#4 plot results####
  
  #hypothesis 1&2: threat and use of arms#### 
  #plot coefficient estimates for change in perceptions based on threat
  
  plot <- matrix(NA,6,3)
  
  plot[,1] <- c(ate.1 - 1.96*se.ate.1, ate.2 - 1.96*se.ate.2, ate.3 - 1.96*se.ate.3,ate.10 - 1.96*se.ate.10, ate.11 - 1.96*se.ate.11, ate.12 - 1.96*se.ate.12)
  plot[,2] <- c(ate.1,ate.2,ate.3,ate.10,ate.11,ate.12)
  plot[,3] <- c(ate.1 + 1.96*se.ate.1, ate.2 + 1.96*se.ate.2, ate.3 + 1.96*se.ate.3,ate.10 + 1.96*se.ate.10, ate.11 + 1.96*se.ate.11, ate.12 + 1.96*se.ate.12)
  
  
  plot<- as.data.frame(plot)
  plot$names <- c("Threat","Threat (Armed Only)","Threat (Unarmed Only)", "Arms","Arms (Threat Only)","Arms (No Threat Only)")
  plot$names <- as.factor(plot$names)
  plot$names <- factor(plot$names, levels = plot$names[7 - row(as.matrix(plot$names))])
  colnames(plot) <- c("lower","ptest","upper","names")
  
  cfs1 <- ggplot(data = plot, aes(x = ptest, y = names, color = names)) +
    geom_point(size=2.5) +
    geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
    geom_vline(xintercept =0, size=1, linetype = "dashed") +
    xlim(-0.25,0.75) + 
    xlab("Change in Proportion Violence") +
    ylab(" ") +
    theme_bw() +
    theme(legend.position = "none") +
    scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50","gray50"))
  
  cfs1
  
  plot2 <- matrix(NA,6,3)
  
  plot2[,1] <- c(ate.4 - 1.96*se.ate.4, ate.5 - 1.96*se.ate.5, ate.6 - 1.96*se.ate.6, ate.13 - 1.96*se.ate.13, ate.14 - 1.96*se.ate.14, ate.15 - 1.96*ate.15)
  plot2[,2] <- c(ate.4,ate.5,ate.6,ate.13,ate.14,ate.15)
  plot2[,3] <- c(ate.4 + 1.96*se.ate.4, ate.5 + 1.96*se.ate.5, ate.6 + 1.96*se.ate.6, ate.13 + 1.96*se.ate.13, ate.14 + 1.96*se.ate.14, ate.15 + 1.96*ate.15)
  
  plot2<- as.data.frame(plot2)
  plot2$names <-  c("Threat","Threat (Armed Only)","Threat (Unarmed Only)", "Arms","Arms (Threat Only)","Arms (No Threat Only)")
  plot2$names <- as.factor(plot2$names)
  plot2$names <- factor(plot2$names, levels = plot2$names[7 - row(as.matrix(plot2$names))])
  colnames(plot2) <- c("lower","ptest","upper","names")
  
  cfs2 <- ggplot(data = plot2, aes(x = ptest, y = names, color = names)) +
    geom_point(size=2.5) +
    geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
    geom_vline(xintercept =0, size=1, linetype = "dashed") +
    xlim(-0.25,0.75) + 
    xlab("Change in Repression Support") +
    ylab(" ") +
    theme_bw() +
    theme(legend.position = "none") +
    scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50","gray50"))
  
  cfs2
  
  #hypothesis 3: in-out group#### 
  #include partisanship interactions
  
  plot4 <- matrix(NA,4,3)
  
  plot4[,1] <- c(model.19int$coefficients[2,1] - 1.96*model.19int$coefficients[2,2], model.19int$coefficients[4,1] - 1.96*model.19int$coefficients[4,2], 
                 model.20int$coefficients[2,1] - 1.96*model.20int$coefficients[2,2], model.20int$coefficients[4,1] - 1.96*model.20int$coefficients[4,2])
  plot4[,2] <- c(model.19int$coefficients[2,1], model.19int$coefficients[4,1],model.20int$coefficients[2,1], model.20int$coefficients[4,1]) 
  plot4[,3] <-  c(model.19int$coefficients[2,1] + 1.96*model.19int$coefficients[2,2], model.19int$coefficients[4,1] + 1.96*model.19int$coefficients[4,2], 
                  model.20int$coefficients[2,1] + 1.96*model.20int$coefficients[2,2], model.20int$coefficients[4,1] + 1.96*model.20int$coefficients[4,2])
  
  plot4<- as.data.frame(plot4)
  plot4$names <- c("BLM vs. Confederate","BLM vs. Con * GOP","BLM vs. All", "BLM vs. All * GOP")
  plot4$names <- as.factor(plot4$names)
  plot4$names <- factor(plot4$names, levels = plot4$names[5 - row(as.matrix(plot4$names))])
  colnames(plot4) <- c("lower","ptest","upper","names")
  
  cfs4 <- ggplot(data = plot4, aes(x = ptest, y = names, color = names)) +
    geom_point(size=2.5) +
    geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
    geom_vline(xintercept =0, size=1, linetype = "dashed") +
    xlim(-0.5,0.75) + 
    xlab("Group Effects on Proportion Violence") +
    ylab(" ") +
    theme_bw() +
    theme(legend.position = "none") +
    scale_color_manual(values=c("gray50","gray50","gray50","gray50"))
 cfs4
  
  plot5 <- matrix(NA,4,3)
  
  plot5[,1] <- c(model.21int$coefficients[2,1] - 1.96*model.21int$coefficients[2,2], model.21int$coefficients[4,1] - 1.96*model.21int$coefficients[4,2], 
                 model.22int$coefficients[2,1] - 1.96*model.22int$coefficients[2,2], model.22int$coefficients[4,1] - 1.96*model.22int$coefficients[4,2])
  plot5[,2] <- c(model.21int$coefficients[2,1], model.21int$coefficients[4,1],model.22int$coefficients[2,1], model.22int$coefficients[4,1]) 
  plot5[,3] <-  c(model.21int$coefficients[2,1] + 1.96*model.21int$coefficients[2,2], model.21int$coefficients[4,1] + 1.96*model.21int$coefficients[4,2], 
                  model.22int$coefficients[2,1] + 1.96*model.22int$coefficients[2,2], model.22int$coefficients[4,1] + 1.96*model.22int$coefficients[4,2])
  
  plot5<- as.data.frame(plot5)
  plot5$names <- c("BLM vs. Confederate","BLM vs. Con * GOP","BLM vs. All", "BLM vs. All * GOP")
  plot5$names <- as.factor(plot5$names)
  plot5$names <- factor(plot5$names, levels = plot5$names[5 - row(as.matrix(plot5$names))])
  colnames(plot5) <- c("lower","ptest","upper","names")
  
  cfs5 <- ggplot(data = plot5, aes(x = ptest, y = names, color = names)) +
    geom_point(size=2.5) +
    geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
    geom_vline(xintercept =0, size=1, linetype = "dashed") +
    xlim(-0.5,0.75) + 
    xlab("Group Effects on Repression Support") +
    ylab(" ") +
    theme_bw() +
    theme(legend.position = "none") +
    scale_color_manual(values=c("gray50","gray50","gray50","gray50")) 
  
  
  cfs5
  
  library(gridExtra)
  
  
  ###################################
  #### Figures 3-4 in Main Text
  ####################################  
  grid.arrange(cfs1,cfs2,ncol=2)
  grid.arrange(cfs4,cfs5,ncol=2)
  
  ############################################################################################################################
  
  
  
#5 tables####
  
  library(stargazer)
  
  #reformat for stargazer
  
  #######################
  ### Table 13 in APpendix
  ####################
  
  model.1 <- lm(action ~ threat, data=samp)
  model.2 <- lm(action ~ threat, data=samp.armedonly)
  model.3 <- lm(action ~ threat, data=samp.unarmedonly)
  
  model.1b <- lm(action ~ threat + age + gender + income + education + pid7, data=samp)
  model.2b <- lm(action ~ threat + age + gender + income + education + pid7, data=samp.armedonly)
  model.3b <- lm(action ~ threat + age + gender + income + education + pid7, data=samp.unarmedonly)
  
  stargazer(model.1, model.2, model.3, model.1b, model.2b, model.3b,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Perception of Violence"),
            covariate.labels=c("Threat of Harm","Age","Gender (Male)", "Gender (Other)","Income","Education","GOP Partisanship"), 
            out="exp_models_1.txt")
  
  #######################
  ### Table 14 in APpendix
  ####################
  
  
  model.4 <- lm(policy ~ threat, data=samp)
  model.5 <- lm(policy ~ threat, data=samp.armedonly)
  model.6 <- lm(policy ~ threat, data=samp.unarmedonly)
  
  model.4b <- lm(policy ~ threat + age + gender + income + education + pid7, data=samp)
  model.5b <- lm(policy ~ threat + age + gender + income + education + pid7, data=samp.armedonly)
  model.6b <- lm(policy ~ threat + age + gender + income + education + pid7, data=samp.unarmedonly)
  
  stargazer(model.4, model.5, model.6, model.4b, model.5b, model.6b,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Support for Repression"),
            covariate.labels=c("Threat of Harm","Age","Gender (Male)", "Gender (Other)","Income","Education","GOP Partisanship"), 
            out="exp_models_2.txt")
  
   
  #######################
  ### Table 15 in APpendix
  ####################
  
  
  model.10 <- lm(action ~ arm, data=samp)
  model.11 <- lm(action ~ arm, data=samp.threatonly)
  model.12 <- lm(action ~ arm, data=samp.unthreatonly)
  
  model.10b <- lm(action ~ arm + age + gender + income + education + pid7, data=samp)
  model.11b <- lm(action ~ arm + age + gender + income + education + pid7, data=samp.threatonly)
  model.12b <- lm(action ~ arm + age + gender + income + education + pid7, data=samp.unthreatonly)
  
  stargazer(model.10, model.11, model.12, model.10b, model.11b, model.12b,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Perception of Violence"),
            covariate.labels=c("Use of Arms","Age","Gender (Male)", "Gender (Other)","Income","Education","GOP Partisanship"), 
            out="exp_models_4.txt")
  
  #######################
  ### Table 16 in APpendix
  ####################
  
  model.13 <- lm(policy ~ arm, data=samp)
  model.14 <- lm(policy ~ arm, data=samp.threatonly)
  model.15 <- lm(policy ~ arm, data=samp.unthreatonly)
  
  model.13b <- lm(policy ~ arm + age + gender + income + education + pid7, data=samp)
  model.14b <- lm(policy ~ arm + age + gender + income + education + pid7, data=samp.threatonly)
  model.15b <- lm(policy ~ arm + age + gender + income + education + pid7, data=samp.unthreatonly)
  
  stargazer(model.13, model.14, model.15, model.13b, model.14b, model.15b,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Support for Repression"),
            covariate.labels=c("Use of Arms","Age","Gender (Male)", "Gender (Other)","Income","Education","GOP Partisanship"), 
            out="exp_models_5.txt")
  
  
  #######################
  ### Table 17 in APpendix
  ####################
  
 model.19int <- lm(action ~ outgrp*gop_id, data=samp.inoutonly)
  model.20int <- lm(action ~ outgrp*gop_id, data=samp)
  model.19intb <- lm(action ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp.inoutonly)
  model.20intb <- lm(action ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp)
  
  stargazer(model.19int, model.20int, model.19intb, model.20intb,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Perception of Violence"),
            covariate.labels=c("Black Lives Matter","GOP (Dummy)","Age","Gender (Male)", "Gender (Other)","Income","Education","GOP Partisanship","BLM*GOP"), 
            out="exp_models_7.txt")
  
  #######################
  ### Table 18 in APpendix
  ####################
  
  
  model.21int <- lm(policy ~ outgrp*gop_id, data=samp.inoutonly)
  model.22int <- lm(policy ~ outgrp*gop_id, data=samp)
  model.21intb <- lm(policy ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp.inoutonly)
  model.22intb <- lm(policy ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp)
  
  stargazer(model.21int, model.22int, model.21intb, model.22intb,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Support for Repression"),
            covariate.labels=c("Black Lives Matter","GOP (Dummy)","Age","Gender (Male)", "Gender (Other)","Income","Education","GOP Partisanship","BLM*GOP"), 
            out="exp_models_8.txt")
  
  #6 testing the mechanism: motive####
  
  #######################
  ### Table 19 in APpendix
  ####################
  
  model.43 <- lm(crimehatemotive ~ threat, data=samp)
  model.44 <- lm(crimehatemotive ~ threat, data=samp.armedonly)
  model.45 <- lm(crimehatemotive ~ threat, data=samp.unarmedonly)
  
  model.43b <- lm(crimehatemotive ~ threat + age + gender + income + education + pid7, data=samp)
  model.44b <- lm(crimehatemotive ~ threat + age + gender + income + education + pid7, data=samp.armedonly)
  model.45b <- lm(crimehatemotive ~ threat + age + gender + income + education + pid7, data=samp.unarmedonly)
  
  stargazer(model.43, model.44, model.45, model.43b, model.44b, model.45b,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Assign as Crime or Hate"),
            covariate.labels=c("Threat of Harm","Age","Gender (Male)", "Gender (Other)","Income","Education","Party ID"), 
            out="exp_models_10.txt")
  
  #######################
  ### Table 20 in APpendix
  ####################
  
  model.46 <- lm(crimehatemotive ~ arm, data=samp)
  model.47 <- lm(crimehatemotive ~ arm, data=samp.threatonly)
  model.48 <- lm(crimehatemotive ~ arm, data=samp.unthreatonly)
  
  model.46b <- lm(crimehatemotive ~ arm + age + gender + income + education + pid7, data=samp)
  model.47b <- lm(crimehatemotive ~ arm + age + gender + income + education + pid7, data=samp.threatonly)
  model.48b <- lm(crimehatemotive ~ arm + age + gender + income + education + pid7, data=samp.unthreatonly)
  
  stargazer(model.46, model.47, model.48, model.46b, model.47b, model.48b,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Assign as Crime or Hate"),
            covariate.labels=c("Use of Arms","Age","Gender (Male)", "Gender (Other)","Income","Education","Party ID"), 
            out="exp_models_11.txt")
  
  #######################
  ### Table 21 in APpendix
  ####################
  
  
  model.49 <- lm(crimehatemotive ~ outgrp*gop_id, data=samp)
  model.50 <- lm(crimehatemotive ~ outgrp*gop_id, data=samp.inoutonly)
  
  model.49b <- lm(crimehatemotive ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp)
  model.50b <- lm(crimehatemotive ~ outgrp*gop_id + age + gender + income + education + pid7, data=samp.inoutonly)
  
  stargazer(model.49, model.50, model.49b, model.50b,
            type="latex", 
            column.separate=NULL,
            dep.var.labels=c("Assign as Crime or Hate"),
            covariate.labels=c("Black Lives Matter","GOP (Dummy)","Age","Gender (Male)", "Gender (Other)","Income","Education","GOP Partisanship","BLM*GOP"), 
            out="exp_models_12.txt")
  ##################################################
  #balance tests, in Appendix, FIgure 2 ####
  #################################################
  
  names <- c("threat","arm","outgrp","education","income","pid7","ideology","age")
  samp.bal <- samp[names]
  samp.bal <- na.omit(samp.bal)
  
  #threat treatment
  
  bal <- matrix(NA,5,3)
  
  bal[,2] <- c((mean(samp.bal$education[samp.bal$threat==1],na.rm=T)-mean(samp.bal$education[samp.bal$threat==0],na.rm=T))/sd(samp.bal$education[samp.bal$threat==1],na.rm=T),
               (mean(samp.bal$income[samp.bal$threat==1],na.rm=T)-mean(samp.bal$income[samp.bal$threat==0],na.rm=T))/sd(samp.bal$income[samp.bal$threat==1],na.rm=T),
               (mean(samp.bal$pid7[samp.bal$threat==1],na.rm=T)-mean(samp.bal$pid7[samp.bal$threat==0],na.rm=T))/sd(samp.bal$pid7[samp.bal$threat==1],na.rm=T),         
               (mean(samp.bal$ideology[samp.bal$threat==1],na.rm=T)-mean(samp.bal$ideology[samp.bal$threat==0],na.rm=T))/sd(samp.bal$ideology[samp.bal$threat==1],na.rm=T),
               (mean(samp.bal$age[samp.bal$threat==1],na.rm=T)-mean(samp.bal$age[samp.bal$threat==0],na.rm=T))/sd(samp.bal$age[samp.bal$threat==1],na.rm=T))
  
  treatnum <- sum(samp.bal$threat==1)
  controlnum <- sum(samp.bal$threat==0)
  
  bal[,1] <- c(bal[1,2] - 1.96*sqrt(var(samp.bal$education[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$threat==1],na.rm=T),
               bal[2,2] - 1.96*sqrt(var(samp.bal$income[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$threat==1],na.rm=T),
               bal[3,2] - 1.96*sqrt(var(samp.bal$pid7[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$pid7[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$pid7[samp.bal$threat==1],na.rm=T),
               bal[4,2] - 1.96*sqrt(var(samp.bal$ideology[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$threat==1],na.rm=T),
               bal[5,2] - 1.96*sqrt(var(samp.bal$age[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$threat==1],na.rm=T))
  
  
  bal[,3] <- c(bal[1,2] + 1.96*sqrt(var(samp.bal$education[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$threat==1],na.rm=T),
               bal[2,2] + 1.96*sqrt(var(samp.bal$income[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$threat==1],na.rm=T),
               bal[3,2] + 1.96*sqrt(var(samp.bal$pid7[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$pid7[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$pid7[samp.bal$threat==1],na.rm=T),
               bal[4,2] + 1.96*sqrt(var(samp.bal$ideology[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$threat==1],na.rm=T),
               bal[5,2] + 1.96*sqrt(var(samp.bal$age[samp.bal$threat==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$threat==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$threat==1],na.rm=T))
  
  bal <- as.data.frame(bal)
  bal$names <- c("Education","Income","Party ID","Ideology","Age")
  colnames(bal) <- c("lower","diff","upper","names")
  
  balplot <- ggplot(data = bal, aes(x = diff, y = names, color = names)) +
    geom_point(size=1.25) +
    geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
    geom_vline(xintercept =0, size=1) +
    xlim(-0.4,0.4) + 
    xlab("Threat Treatments") +
    ylab("Covariates") +
    theme_bw() + 
    theme(legend.position = "none") +
    scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50"))
  balplot
  #arms treatment
  
  bal2 <- matrix(NA,5,3)
  
  bal2[,2] <- c((mean(samp.bal$education[samp.bal$arm==1],na.rm=T)-mean(samp.bal$education[samp.bal$arm==0],na.rm=T))/sd(samp.bal$education[samp.bal$arm==1],na.rm=T),
                (mean(samp.bal$income[samp.bal$arm==1],na.rm=T)-mean(samp.bal$income[samp.bal$arm==0],na.rm=T))/sd(samp.bal$income[samp.bal$arm==1],na.rm=T),
                (mean(samp.bal$pid7[samp.bal$arm==1],na.rm=T)-mean(samp.bal$pid7[samp.bal$arm==0],na.rm=T))/sd(samp.bal$pid7[samp.bal$arm==1],na.rm=T),         
                (mean(samp.bal$ideology[samp.bal$arm==1],na.rm=T)-mean(samp.bal$ideology[samp.bal$arm==0],na.rm=T))/sd(samp.bal$ideology[samp.bal$arm==1],na.rm=T),
                (mean(samp.bal$age[samp.bal$arm==1],na.rm=T)-mean(samp.bal$age[samp.bal$arm==0],na.rm=T))/sd(samp.bal$age[samp.bal$arm==1],na.rm=T))
  
  treatnum <- sum(samp.bal$arm==1)
  controlnum <- sum(samp.bal$arm==0)
  
  bal2[,1] <- c(bal2[1,2] - 1.96*sqrt(var(samp.bal$education[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$arm==1],na.rm=T),
                bal2[2,2] - 1.96*sqrt(var(samp.bal$income[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$arm==1],na.rm=T),
                bal2[3,2] - 1.96*sqrt(var(samp.bal$pid7[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$pid7[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$pid7[samp.bal$arm==1],na.rm=T),
                bal2[4,2] - 1.96*sqrt(var(samp.bal$ideology[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$arm==1],na.rm=T),
                bal2[5,2] - 1.96*sqrt(var(samp.bal$age[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$arm==1],na.rm=T))
  
  
  bal2[,3] <- c(bal2[1,2] + 1.96*sqrt(var(samp.bal$education[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$arm==1],na.rm=T),
                bal2[2,2] + 1.96*sqrt(var(samp.bal$income[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$arm==1],na.rm=T),
                bal2[3,2] + 1.96*sqrt(var(samp.bal$pid7[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$pid7[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$pid7[samp.bal$arm==1],na.rm=T),
                bal2[4,2] + 1.96*sqrt(var(samp.bal$ideology[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$arm==1],na.rm=T),
                bal2[5,2] + 1.96*sqrt(var(samp.bal$age[samp.bal$arm==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$arm==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$arm==1],na.rm=T))
  
  bal2 <- as.data.frame(bal2)
  bal2$names <- c("Education","Income","Party ID","Ideology","Age")
  colnames(bal2) <- c("lower","diff","upper","names")
  
  balplot2 <- ggplot(data = bal2, aes(x = diff, y = names, color = names)) +
    geom_point(size=1.25) +
    geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
    geom_vline(xintercept =0, size=1) +
    xlim(-0.4,0.4) + 
    xlab("Arms Treatment") +
    ylab(" ") +
    theme_bw() +
    theme(legend.position = "none")+
    scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50"))
  balplot2
  
  
  #outgroup treatment
  
  bal3 <- matrix(NA,5,3)
  
  bal3[,2] <- c((mean(samp.bal$education[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$education[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$education[samp.bal$outgrp==1],na.rm=T),
                (mean(samp.bal$income[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$income[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$income[samp.bal$outgrp==1],na.rm=T),
                (mean(samp.bal$pid7[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$pid7[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$pid7[samp.bal$outgrp==1],na.rm=T),         
                (mean(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$ideology[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T),
                (mean(samp.bal$age[samp.bal$outgrp==1],na.rm=T)-mean(samp.bal$age[samp.bal$outgrp==0],na.rm=T))/sd(samp.bal$age[samp.bal$outgrp==1],na.rm=T))
  
  treatnum <- sum(samp.bal$outgrp==1)
  controlnum <- sum(samp.bal$outgrp==0)
  
  bal3[,1] <- c(bal3[1,2] - 1.96*sqrt(var(samp.bal$education[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$outgrp==1],na.rm=T),
                bal3[2,2] - 1.96*sqrt(var(samp.bal$income[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$outgrp==1],na.rm=T),
                bal3[3,2] - 1.96*sqrt(var(samp.bal$pid7[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$pid7[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$pid7[samp.bal$outgrp==1],na.rm=T),
                bal3[4,2] - 1.96*sqrt(var(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T),
                bal3[5,2] - 1.96*sqrt(var(samp.bal$age[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$outgrp==1],na.rm=T))
  
  
  bal3[,3] <- c(bal3[1,2] + 1.96*sqrt(var(samp.bal$education[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$education[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$education[samp.bal$outgrp==1],na.rm=T),
                bal3[2,2] + 1.96*sqrt(var(samp.bal$income[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$income[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$income[samp.bal$outgrp==1],na.rm=T),
                bal3[3,2] + 1.96*sqrt(var(samp.bal$pid7[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$pid7[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$pid7[samp.bal$outgrp==1],na.rm=T),
                bal3[4,2] + 1.96*sqrt(var(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$ideology[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$ideology[samp.bal$outgrp==1],na.rm=T),
                bal3[5,2] + 1.96*sqrt(var(samp.bal$age[samp.bal$outgrp==1],na.rm=T)/treatnum + var(samp.bal$age[samp.bal$outgrp==0],na.rm=T)/controlnum)/sd(samp.bal$age[samp.bal$outgrp==1],na.rm=T))
  
  bal3 <- as.data.frame(bal3)
  bal3$names <- c("Education","Income","Party ID","Ideology","Age")
  colnames(bal3) <- c("lower","diff","upper","names")
  
  balplot3 <- ggplot(data = bal3, aes(x = diff, y = names, color = names)) +
    geom_point(size=1.25) +
    geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
    geom_vline(xintercept =0, size=1) +
    xlim(-0.4,0.4) + 
    xlab("Out-Group Treatment") +
    ylab(" ") +
    theme_bw() +
    theme(legend.position = "none") +
    scale_color_manual(values=c("gray50","gray50","gray50","gray50","gray50"))
  
  library(gridExtra)
  
  ###################################
  #### Figure 2 in Appendix
  ####################################
  grid.arrange(balplot,balplot2,balplot3, ncol=3)
  
  #############################################
  #attrition check####
  ##############################################
  #MIPO?
  attrition <- summary(lm(manipulation ~ factor(treatment) + gender + age + education + income + pid7, data=samp))
  
  #test with randomizataion inference according to Gerber and Green
  ri.lm.multiple <- function(df, outcome, covariate1, covariate2, covariate3, covariate4, covariate5, reps, n, m, num){
    
    library(randomizr)
    #create objects
    Fs <- vector(mode="numeric")
    assns <- as.data.frame(matrix(NA,nrow=nrow(df),ncol=reps))
    dfs <- list()
    
    #generate random assignments
    for(i in 1:reps){
      assns[,i] <- complete_ra(N = n, prob_each = m, num_arms = num)}
    
    #remove duplicate assignments
    n=reps
    assns <- assns[!duplicated(lapply(assns, function(k) head(k, n)))]
    
    #take F statistic from regression and store
    for(i in 1:ncol(assns)){dfs[[i]] <- cbind(outcome, assns[,i], covariate1,covariate2, covariate3, covariate4, covariate5) 
    Fs[i] <- summary(lm(dfs[[i]][,1] ~ factor(dfs[[i]][,2]) + dfs[[i]][,3]+ dfs[[i]][,4]+ dfs[[i]][,5]+ 
                          dfs[[i]][,6]+ dfs[[i]][,7]))$fstatistic[1]}
    Fs <- as.data.frame(sort(Fs))
    
    #output
    return(Fs)
  }
  
  fstats <- ri.lm.multiple(samp, samp$manipulation, samp$gender, samp$age, samp$education, samp$income, samp$pid7,
                           1000,1137,rep(1/12,12),12)
  
  sig <- sum(ifelse(fstats >= attrition$fstatistic[1],1,0))
  
  #####################################################################################
  ############END OF US SURVEY
  ############# MAKE SURE TO CLEAR VARIABLE NAMES BETWEEN US AND ISRAEL SURVEY
  ############# CODE WILL BREAK OTHERWISE, DUE TO NAME REPITIONS
  #################################################################################
  
  #################################
  # Data Comparisons - Figure 1 in Appendix 
  ################################
  rm(list = ls())
  cat("\014")
  set.seed(30030)
  
  #1 load data (us)####
  compare <- read.csv("US_Israel_compare.csv")
  #4 plot descriptive comparisons####
  
  library(ggplot2)
  library(scales)
  library(gridExtra)
  
  hists1 <- ggplot(compare,aes(ideology,..density..,colour=group)) +
    geom_freqpoly(binwidth=1, size=1.5) +
    coord_cartesian(xlim=c(1,9)) +
    xlab("Ideology") +
    ylab("Density") +
    scale_color_manual(values=c('steelblue3','firebrick3'))+
    theme(panel.background = element_blank(), 
          text = element_text(size=14),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          legend.title=element_blank())
  
  hists2 <- ggplot(compare,aes(education,..density..,colour=group)) +
    geom_freqpoly(binwidth=1, size=1.5) +
    coord_cartesian(xlim=c(1,5)) +
    xlab("Education") +
    ylab("Density") +
    scale_color_manual(values=c('steelblue3','firebrick3'))+
    theme(panel.background = element_blank(), 
          text = element_text(size=14),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          legend.title=element_blank())
  
  hists3 <- ggplot(compare,aes(income,..density..,colour=group)) +
    geom_freqpoly(binwidth=1, size=1.5) +
    coord_cartesian(xlim=c(1,5)) +
    xlab("Income") +
    ylab("Density") +
    scale_color_manual(values=c('steelblue3','firebrick3'))+
    theme(panel.background = element_blank(), 
          text = element_text(size=14),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          legend.title=element_blank())
  
  hists4 <- ggplot(compare,aes(age,..density..,colour=group)) +
    geom_freqpoly(binwidth=5, size=1.5) +
    coord_cartesian(xlim=c(18,90)) +
    xlab("Age") +
    ylab("Density") +
    scale_color_manual(values=c('steelblue3','firebrick3'))+
    theme(panel.background = element_blank(), 
          text = element_text(size=14),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          legend.title=element_blank())
  
  compare <- subset(compare, !is.na(compare$newssrc))
  compare <- subset(compare, !is.na(compare$gender))
  
  hists5 <- ggplot(compare,aes(newssrc,group=group)) +
    geom_bar(aes(y = ..prop.., fill = group), position="dodge", stat="count") + 
    scale_y_continuous(labels=scales::percent) +
    xlab("News Source") +
    ylab("Percent") +
    scale_fill_manual(values=c('steelblue3','firebrick3'))+
    theme(panel.background = element_blank(), 
          text = element_text(size=14),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          legend.title=element_blank())
  
  hists6 <- ggplot(compare,aes(gender,group=group)) +
    geom_bar(aes(y = ..prop.., fill = group), position="dodge", stat="count") + 
    scale_y_continuous(labels=scales::percent) +
    xlab("Gender") +
    ylab("Percent") +
    scale_fill_manual(values=c('steelblue3','firebrick3'))+
    theme(panel.background = element_blank(), 
          text = element_text(size=14),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          legend.title=element_blank())
  
  grid.arrange(hists1,hists2,hists3,hists4,hists5,hists6, nrow=2, ncol=3)
  